import sys
# Our numerical workhorses
import numpy as np
import pandas as pd
import scipy.optimize
from scipy.integrate import quad
import scipy.stats as st
# Numerical differentiation package
import numdifftools as ndt
# MCMC package
import emcee
# MCMC results visualization package
import corner
# Import pyplot for plotting
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LightSource
# Seaborn, useful for graphics
#import seaborn as sns
# Bokeh stuff
import bokeh.charts
import bokeh.io
# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline
# This enables high res graphics inline (only use with static plots (non-Bokeh))
# SVG is preferred, but there is a bug in Jupyter with vertical lines
%config InlineBackend.figure_formats = {'svg', 'retina'}
# JB's favorite Seaborn settings for notebooks
rc = {'lines.linewidth': 2,
'axes.labelsize': 18,
'axes.titlesize': 18,
'axes.facecolor': 'DFDFE5'}
#sns.set_context('notebook', font_scale=1.5, rc={'lines.linewidth': 2.5})
#sns.set_style('darkgrid', {'axes.facecolor': '(0.875, 0.875, 0.9)'})
# Display graphics in this notebook
bokeh.io.output_notebook()
df = pd.read_excel('data/Copy of reversible lambda for Jihoon.xlsx', comment='#', )
df
df['l1'] = df['l1'] + 50.0
df['l2'] = df['l2'] + 50.0
const_drug = df[:9]
pulse_drug = df.drop(df.index[6:10])
pulse_drug = pulse_drug.drop(pulse_drug.index[0:5])
const_drug = const_drug.sort('days after starting drug treatment')
pulse_drug = pulse_drug.sort('days after starting drug treatment')
const_drug_alt = const_drug
const_drug_alt
# Add a few dummy data points at the end to try to induce a steady state equilibrium and deter oscillation
# to pulsed case
pulse_drug = pulse_drug.append(pd.DataFrame([['drug off', 'dr29.D40', 69, 88.022285, 12.981519], \
['drug off', 'dr29.D45', 74, 88.022285, 12.981519], \
['drug off', 'dr29.D50', 79, 88.022285, 12.981519]], \
columns=['drug condition', 'name', 'days after starting drug treatment', 'l1', 'l2']))
pulse_drug = pulse_drug.reset_index()
pulse_drug = pulse_drug.drop('index', 1)
const_drug = pulse_drug
const_drug
# Theoretical model for populations
def popDynamics(p, d):
data_times, (c01, c02) = d
data_times = np.array(data_times)
a1, a2, a3, alpha1, beta1, alpha2, s = p
euler_step = 0.5
td = np.arange(data_times[0], data_times[-1] + euler_step, euler_step)
l_theor_full = np.zeros([len(td), len(d[1])])
l_theor_full[0] = np.array([[c01, c02]])
l_theor = np.array([[c01, c02]])
# Temporarily? changed the diff eqs such that a1 and a2 are the same, in order to make the resulting vector field
# conservative --> more biologically plausible
for i in range(1, len(td)):
l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
- a1 * l_theor_full[i-1][1])
l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
- a3 * l_theor_full[i-1][1])
if td[i] in data_times:
l_theor = np.append(l_theor, [[l_theor_full[i][0], l_theor_full[i][1]]], axis=0)
return l_theor_full, l_theor
def log_posterior(p, d, l):
'''
Assuming uniform priors for all of the parameters and Jeffreys prior for sigma. Assuming Gaussian
distributed measurements in the experiment.
'''
data_times, (c01, c02) = d
a1, a2, a3, alpha1, beta1, alpha2, s = p
# Zero probability of having non-positive value of stdev
if s <= 0:
return -np.inf
_, pops_theor = popDynamics(p, d)
# Zero probability of having absurd population values
for i in range(len(pops_theor)):
for j in range(2):
if np.abs(pops_theor[i][j]) > 300.0 or np.isnan(pops_theor[i][j]) or np.abs(pops_theor[i][j]) < 0:
return -np.inf
for val in p:
if val < 0:
return -np.inf
#if alpha1 < 0 or alpha2 < 0 or a2 < 0 or a1 < 0:
#return -np.inf
if a1 > 1 or a2 > 1 or a3 > 1:
return -np.inf
# Nyquist frequency: 1 day^-1 so keep frequency below or equal to that
if 1 / (2 * np.pi) * np.sqrt(np.sqrt(a1**2.0 * a3**2.0) + np.sqrt(beta1**2.0 * a2**2.0)) > 1.0:
return -np.inf
return -(2 * len(data_times) + 1) * np.log(s) - 0.5 / s**2.0 * (np.sum((l[:,1]-pops_theor[:,1])**2.0) \
+ np.sum((l[:,0]-pops_theor[:,0])**2.0))
n_dim = 7 # 7 parameters in the model
n_walkers = 1000 # Number of MCMC walkers
n_burn = 500 # 500 burn-in steps
n_steps = 5000 # Total number of steps after burn-in
# Seeding random number generator
np.random.seed(1)
# p0[i, j] is the starting point for the ith walk for the jth variable
p0 = np.empty((n_walkers, n_dim))
for i in range(n_dim - 1):
p0[:, i] = np.random.uniform(0, 1, n_walkers)
p0[:, n_dim-1] = np.random.exponential(0.1, n_walkers) #sigma
# Creating a MCMC sampler
emp_data = np.array([const_drug['l1'], const_drug['l2']])
emp_data = np.transpose(emp_data)
sampler = emcee.EnsembleSampler(n_walkers, n_dim, log_posterior, args=((const_drug['days after starting drug treatment'], \
emp_data[0]), emp_data,), threads=1)
# Using 1 core since multithreading seems to be broken
# Do burn-in
pos, prob, state = sampler.run_mcmc(p0, n_burn, storechain=False)
#Actual run
for i, result in enumerate(sampler.sample(pos, lnprob0=prob, iterations=n_steps)):
if (i+1) % 100 == 0:
print("{0:5.1%}".format(float(i) / n_steps))
fig, ax = plt.subplots(7, 1, sharex=True)
for i in range(6):
ax[i].plot(sampler.chain[0,:,i], 'k-', lw=0.4)
ax[i].plot([0, n_steps-1],
[sampler.chain[0,:,i].mean(), sampler.chain[0,:,i].mean()], 'r-', lw=0.4)
ax[6].set_xlabel('sample number')
# Get the index of the most probable parameter set
max_ind = np.argmax(sampler.flatlnprobability)
# Pull out values.
a1, a2, a3, alpha1, beta1, alpha2, s = sampler.flatchain[max_ind,:]
# Calculate errors
a1_e, a2_e, a3_e, alpha1_e, beta1_e, alpha2_e, s_e = sampler.flatchain.std(axis=0)
# Print the results
print('Most probable parameter values:' + \
'\na1: ' + str(a1) + ' +- ' + str(a1_e) + \
'\na2: ' + str(a2) + ' +- ' + str(a2_e) + \
'\na3: ' + str(a3) + ' +- ' + str(a3_e) + \
'\nalpha1: ' + str(alpha1) + ' +- ' + str(alpha1_e) + \
'\nbeta1: ' + str(beta1) + ' +- ' + str(beta1_e) + \
'\nalpha2: ' + str(alpha2) + ' +- ' + str(alpha2_e) + \
'\ns: ' + str(s) + ' +- ' + str(s_e))
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, [0.0, 0.0]))[0], alpha=0.8)
_ = plt.scatter(const_drug['days after starting drug treatment'], const_drug['l2'], color='green')
_ = plt.scatter(const_drug['days after starting drug treatment'], const_drug['l1'], color='blue')
euler_step = 0.5
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0], lw=0.6, alpha=0.8)
_ = plt.plot(const_drug['days after starting drug treatment'], \
popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (const_drug['days after starting drug treatment'], emp_data[0]))[1], lw=1.5)
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
np.array(const_drug['days after starting drug treatment'])[-1] + 100 * euler_step, euler_step)
_ = plt.plot(td, popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0])
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
np.savetxt('data/pulse_drug_results.csv', popDynamics((a1, a2, a3, alpha1, beta1, alpha2, s), (td, emp_data[0]))[0], \
delimiter=',')
# Plotting the probability distributions of the most probable parameter values
#fig = corner.corner(sampler.flatchain, labels=[r'$a_1$', r'$a_2$', r'$a_3$', r'$\alpha_1$', r'$\beta_1$', \
# r'$\alpha_2$', r'$\sigma$'], bins=100, lw=2)
#fig.savefig("data/corner_pulse.png")
Now we do the same for the constant drug (no release):
euler_step = 0.5
emp_data_alt = np.array([const_drug_alt['l1'], const_drug_alt['l2']])
emp_data_alt = np.transpose(emp_data_alt)
sampler_const = emcee.EnsembleSampler(n_walkers, n_dim, log_posterior, args=((const_drug_alt['days after starting drug treatment'], \
emp_data_alt[0]), emp_data_alt,), threads=1)
# Using 1 core since multithreading seems to be broken
# Do burn-in
pos, prob, state = sampler_const.run_mcmc(p0, n_burn, storechain=False)
for i, result in enumerate(sampler_const.sample(pos, lnprob0=prob, iterations=n_steps)):
if (i+1) % 100 == 0:
print("{0:5.1%}".format(float(i) / n_steps))
# Get the index of the most probable parameter set
max_ind = np.argmax(sampler_const.flatlnprobability)
# Pull out values.
a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc = sampler_const.flatchain[max_ind,:]
# Calculate errors
a1c_e, a2c_e, a3c_e, alpha1c_e, beta1c_e, alpha2c_e, sc_e = sampler_const.flatchain.std(axis=0)
# Print the results
print('Most probable parameter values:' + \
'\na1: ' + str(a1c) + ' +- ' + str(a1c_e) + \
'\na2: ' + str(a2c) + ' +- ' + str(a2c_e) + \
'\na3: ' + str(a3c) + ' +- ' + str(a3c_e) + \
'\nalpha1: ' + str(alpha1c) + ' +- ' + str(alpha1c_e) + \
'\nbeta1: ' + str(beta1c) + ' +- ' + str(beta1c_e) + \
'\nalpha2: ' + str(alpha2c) + ' +- ' + str(alpha2c_e) + \
'\ns: ' + str(sc) + ' +- ' + str(sc_e))
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
np.array(const_drug_alt['days after starting drug treatment'])[-1] + 100 * euler_step, euler_step)
_ = plt.plot(td, popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0])
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
np.savetxt('data/const_drug_results.csv', popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0], \
delimiter=',')
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l2'], color='green')
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l1'], color='blue')
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(td, popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (td, emp_data_alt[0]))[0], lw=0.6, alpha=0.8)
_ = plt.plot(const_drug_alt['days after starting drug treatment'], \
popDynamics((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), (const_drug_alt['days after starting drug treatment'], emp_data_alt[0]))[1], lw=1.5)
_ = plt.legend(['l1', 'l2'], loc='upper right')
_ = plt.xlabel('Cycle')
_ = plt.ylabel('population')
_ = plt.show()
# Plotting the probability distributions of the most probable parameter values CONSTANT DRUG
#fig = corner.corner(sampler_const.flatchain, labels=[r'$a_1$', r'$a_3$', r'$\alpha_1$', r'$\beta_1$', \
# r'$\alpha_2$', r'$\sigma$'], bins=100, lw=2)
#fig.savefig('data/corner_const.png')
def distances(p, l1Domain, l2Domain, td, pulseTrue):
# Calculating trajectory length from start to steady state
# l1Domain and l2Domain are the starting domains of l1 and l2. Must be equally sized
# Overlays the trajectory using empirical starting points as well, with z-value equal to distance left to cover.
if pulseTrue:
dist_emp_data = emp_data[0]
else:
dist_emp_data = emp_data_alt[0]
a1, a2, a3, alpha1, beta1, alpha2, s = p
euler_step = 0.5
dist_comb = 0
l_theor_full = np.zeros([len(td), 2, len(l1Domain), len(l2Domain)])
l_theor_full[0] = np.meshgrid(l1Domain, l2Domain)
for i in range(1, len(td)):
l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
- a1 * l_theor_full[i-1][1])
l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
- a3 * l_theor_full[i-1][1])
dist_1 = l_theor_full[i][0] - l_theor_full[i-1][0]
dist_2 = l_theor_full[i][1] - l_theor_full[i-1][1]
dist_comb = dist_comb + np.sqrt(dist_1**2.0 + dist_2**2.0)
l1Domain, l2Domain = np.meshgrid(l1Domain, l2Domain)
pops = popDynamics(p, (td, dist_emp_data))[0]
dist_traj_indiv = np.sqrt(np.diff([val[0] for val in pops])**2.0 + np.diff([val[1] for val in pops])**2.0)
dist_traj_indiv = np.append(dist_traj_indiv, dist_traj_indiv[-1])
dist_traj = np.zeros(len(dist_traj_indiv))
for i in range(len(dist_traj_indiv)):
dist_traj[i] = np.sum(dist_traj_indiv) - np.sum(dist_traj_indiv[:i])
#dist_traj = dist_traj[::-1]
UN = alpha1 - beta1 * l1Domain - a1 * l2Domain
VN = alpha2 + a2 * l1Domain - a3 * l2Domain
fig = plt.figure(figsize=(12,9))
ax = fig.gca(projection='3d')
pulse_surf = ax.plot_surface(l1Domain, l2Domain, dist_comb, cstride=10, rstride=10, cmap=cm.jet, \
linewidth=0, antialiased=False, alpha=0.6)
fig.colorbar(pulse_surf, shrink=0.5, aspect=5)
pulse_surf = ax.quiver(l1Domain[::20, ::20], l2Domain[::20, ::20], dist_comb[::20, ::20], UN[::20, ::20], \
VN[::20, ::20], 0, length=1.0, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
pulse_surf = ax.plot([val[0] for val in pops], [val[1] for val in pops], dist_traj, lw=1.4, color='white', \
label='From exp. starting point')
pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
pulse_surf = ax.set_zlabel('Total distance of trajectory to steady state')
if pulseTrue == False:
ax.view_init(azim=200, elev=30)
ax.set_title('Constant drug condition trajectory lengths')
else:
ax.set_title('Pulse drug condition trajectory lengths')
ax.view_init(azim=-40, elev=30)
return np.min(np.min(dist_comb))
s = 2.0
euler_step = 0.5
# Pulse drug
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances((a1, a2, a3, alpha1, beta1, alpha2, s), l1Domain, l2Domain, td, True)
sc = 2.0
# Constant drug
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
distances((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), l1Domain, l2Domain, td, False)
def distances_2d(p, l1Domain, l2Domain, td, pulseTrue):
# Calculating trajectory length from start to steady state
# l1Domain and l2Domain are the starting domains of l1 and l2. Must be equally sized
# Overlays the trajectory using empirical starting points as well
if pulseTrue:
dist_emp_data = emp_data[0]
else:
dist_emp_data = emp_data_alt[0]
a1, a2, a3, alpha1, beta1, alpha2, s = p
euler_step = 0.5
dist_comb = 0
l_theor_full = np.zeros([len(td), 2, len(l1Domain), len(l2Domain)])
l_theor_full[0] = np.meshgrid(l1Domain, l2Domain)
for i in range(1, len(td)):
l_theor_full[i][0] = l_theor_full[i-1][0] + euler_step * (alpha1 - beta1 * l_theor_full[i-1][0] \
- a1 * l_theor_full[i-1][1])
l_theor_full[i][1] = l_theor_full[i-1][1] + euler_step * (alpha2 + a2 * l_theor_full[i-1][0] \
- a3 * l_theor_full[i-1][1])
dist_1 = l_theor_full[i][0] - l_theor_full[i-1][0]
dist_2 = l_theor_full[i][1] - l_theor_full[i-1][1]
dist_comb = dist_comb + np.sqrt(dist_1**2.0 + dist_2**2.0)
l1Domain, l2Domain = np.meshgrid(l1Domain, l2Domain)
pops = popDynamics(p, (td, dist_emp_data))[0]
dist_traj_indiv = np.sqrt(np.diff([val[0] for val in pops])**2.0 + np.diff([val[1] for val in pops])**2.0)
dist_traj_indiv = np.append(dist_traj_indiv, dist_traj_indiv[-1])
dist_traj = np.zeros(len(dist_traj_indiv))
for i in range(len(dist_traj_indiv)):
dist_traj[i] = np.sum(dist_traj_indiv) - np.sum(dist_traj_indiv[:i])
UN = alpha1 - beta1 * l1Domain - a1 * l2Domain
VN = alpha2 + a2 * l1Domain - a3 * l2Domain
fig = plt.figure(figsize=(12,9))
ax = fig.gca()
pulse_surf = plt.plot([val[0] for val in pops], [val[1] for val in pops], lw=1.4, color='white')
pulse_surf = plt.pcolormesh(l1Domain, l2Domain, dist_comb, cmap=cm.jet, \
linewidth=0, antialiased=False, alpha=0.6)
fig.colorbar(pulse_surf, shrink=0.5, aspect=5, label='Trajectory distance')
pulse_surf = plt.quiver(l1Domain[::20, ::20], l2Domain[::20, ::20], UN[::20, ::20], \
VN[::20, ::20], color='black', pivot='middle')
pulse_surf = plt.xlabel(r'$\lambda_1$ starting value')
pulse_surf = plt.ylabel(r'$\lambda_2$ starting value')
if pulseTrue == False:
plt.title('Constant drug condition trajectory lengths')
else:
plt.title('Pulse drug condition trajectory lengths')
return np.min(np.min(dist_comb))
# Pulse drug
td = np.arange(np.array(const_drug['days after starting drug treatment'])[0], \
np.array(const_drug['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances_2d((a1, a2, a3, alpha1, beta1, alpha2, s), l1Domain, l2Domain, td, True)
# Constant drug
# Pulse drug
td = np.arange(np.array(const_drug_alt['days after starting drug treatment'])[0], \
np.array(const_drug_alt['days after starting drug treatment'])[-1] + euler_step, euler_step)
l1Domain = np.arange(-30, 210, 1)
l2Domain = np.arange(-30, 210, 1)
distances_2d((a1c, a2c, a3c, alpha1c, beta1c, alpha2c, sc), l1Domain, l2Domain, td, False)
# More analysis of previously found model fitting
df_pulse = pd.read_csv('data/pulse_drug_results.csv', comment='#', )
df_const = pd.read_csv('data/const_drug_results.csv', comment='#',)
df_pulse.head()
_ = plt.scatter(pulse_drug['days after starting drug treatment'], pulse_drug['l2']-50.0, color='green')
_ = plt.scatter(pulse_drug['days after starting drug treatment'], pulse_drug['l1']-50.0, color='blue')
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(df_pulse['time (cycles)'], df_pulse['l1']-50.0, df_pulse['time (cycles)'], df_pulse['l2']-50.0)
_ = plt.legend([r'$\lambda_1$', r'$\lambda_2$'], loc='upper right')
_ = plt.xlabel('Day')
_ = plt.ylabel('population')
_ = plt.show()
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l2']-50.0, color='green')
_ = plt.scatter(const_drug_alt['days after starting drug treatment'], const_drug_alt['l1']-50.0, color='blue')
plt.gca().set_color_cycle(['blue', 'green', 'blue', 'green'])
_ = plt.plot(df_const['time (cycle)'], df_const['l1']-50.0, df_const['time (cycle)'], df_const['l2']-50.0,)
_ = plt.legend([r'$\lambda_1$', r'$\lambda_2$'], loc='upper right')
_ = plt.xlabel('Day')
_ = plt.ylabel('population')
_ = plt.show()
df_landscape = pd.read_excel('data/input_reversible_landscape.xlsx', comment='#')
df_landscape.head()
xi, yi = np.meshgrid(df_landscape['xi'].unique() + 50.0, df_landscape['yi'].unique() + 50.0)
z = np.zeros((len(xi), len(yi)))
for i in range(len(xi)):
for j in range(len(yi)):
z[i][j] = df_landscape['z'][j*len(xi) + i] / 1000.0
a1 = 0.248
a2 = 0.016
a3 = 0.109
alpha1 = 11.665
beta1 = 0.090
alpha2 = 0.124
UN = alpha1 - beta1 * xi - a1 * yi
VN = alpha2 + a2 * xi - a3 * yi
# Plotting landscape
fig = plt.figure(figsize=(12, 8))
ax = fig.gca(projection='3d')
# Create light source object.
ls = LightSource(azdeg=285, altdeg=70, hsv_min_val=0.6, hsv_min_sat=0.8)
# Shade data, creating an rgb array.
rgb = ls.shade(z, plt.cm.rainbow)
pulse_surf = ax.plot_surface(xi-50, yi-50, z, cstride=1, rstride=1, facecolors=rgb, \
linewidth=0, antialiased=True, alpha=1, shade=True)
pulse_surf = ax.quiver(xi-50, yi-50, z, UN, \
VN, 0, length=0.65, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
pulse_surf = ax.set_zlabel(r'Free energy ($\times10^3$)')
ax.view_init(azim=240, elev=60)
#ax.dist=0.1
plt.title('Pulse drug condition energy landscape')
xi, yi = np.meshgrid(df_landscape['xi'].unique() + 50.0, df_landscape['yi'].unique() + 50.0)
z = np.zeros((len(xi), len(yi)))
for i in range(len(xi)):
for j in range(len(yi)):
z[i][j] = df_landscape['z'][j*len(xi) + i] / 1000.0
a1c = 0.339
a2c = 0.076
a3c = 0.400
alpha1c = 24.181
beta1c = 0.0003
alpha2c = 28.618
UN = alpha1c - beta1c * xi - a1c * yi
VN = alpha2c + a2c * xi - a3c * yi
# Plotting landscape
fig = plt.figure(figsize=(12, 8))
ax = fig.gca(projection='3d')
# Create light source object.
ls = LightSource(azdeg=285, altdeg=70, hsv_min_val=0.6, hsv_min_sat=0.8)
# Shade data, creating an rgb array.
rgb = ls.shade(z, plt.cm.rainbow)
pulse_surf = ax.plot_surface(xi-50, yi-50, z, cstride=1, rstride=1, facecolors=rgb, \
linewidth=0, antialiased=True, alpha=1, shade=True)
pulse_surf = ax.quiver(xi-50, yi-50, z, UN, \
VN, 0, length=0.65, linewidths=1.5, arrow_length_ratio=0.6, color='black', pivot='middle')
pulse_surf = ax.set_xlabel(r'$\lambda_1$ starting value')
pulse_surf = ax.set_ylabel(r'$\lambda_2$ starting value')
pulse_surf = ax.set_zlabel(r'Free energy ($\times10^3$)')
ax.view_init(azim=240, elev=60)
#ax.dist=0.1
plt.title('Constant drug condition energy landscape')